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We explicitly calculate the lowest order systematic inhomogeneity-induced corrections to the cos- 
mological constant that one would infer from an analysis of the luminosities and redshifts of Type 
la supernovae, assuming a homogeneous universe. The calculation entails a post-Newtonian expan- 
sion within the framework of second order perturbation theory, wherein we consider the effects of 
subhorizon density perturbations in a flat, dust dominated universe. Within this formalism, we 
calculate luminosity distances and redshifts along the past light cone of an observer. The luminosity 
distance- redshift relation is then averaged over viewing angles and ensemble averaged, assuming that 
density fluctuations at a given cosmic time are a homogeneous random process. The resulting rela- 
tion is fit to that of a homogeneous model containing dust and a cosmological constant, in order to 
deduce the best-fit cosmological constant density S7a- We find that the luminosity distance-redshift 
relation is indeed modified, even for large sample sizes, but only by a very small fraction, of order 
10 _J for z ~ 0.1. This lowest order deviation depends on the peculiar velocities of the source and 
the observer. However, when fitting this perturbed relation to that of a homogeneous universe, via 
maximizing a likelihood function, we find that the inferred cosmological constant can be surpris- 
ingly large, depending on the range of redshifts sampled. For a sample of supernovae extending from 
Zmin ~ 0.02 out to a limiting redshift z max = 0.15, we find that Q.a ~ 0.004. The value of Qa has 
a large variance, and its magnitude tends to get progressively larger as the limiting redshift Zmax 
gets smaller, implying that precision measurements of Qa from nearby supernova data will require 
taking this effect into account. This effect has been referred to in the past as the "fitting problem", 
and more recently as subhorizon "backreaction" . We find that it is likely too small to explain the 
observed value Qa ~ 0.7. There have been previous claims of much larger backreaction effects. By 
contrast to those calculations, our work is directly related to how observers deduce cosmological 
parameters from astronomical data. 

PACS numbers: 98.80.-k, 98.80.Jk, 98.80.Es 

I. INTRODUCTION 

It appears as though the Universe is expanding at an accelerating rate, as has been deduced from luminosity distance 
measurements of Type la supernovae, which appear dimmer than one would expect based on general relativity without 
a cosmological constant [H, H|. This acceleration has also been deduced from measurements of the current matter 
density Qm rj 0.27, which is too small to close the Universe as required by cosmic microwave background radiation 
(CMB) observations with Hq priors from HST [1,0]. Explanations for this discrepancy have been put forward, but 
most employ a modification of general relativity on cosmological scales or the addition of an exotic "dark energy" 
field. 

There have also been attempts to explain this seemingly anomalous cosmic acceleration as a consequence of sub- 
horizon inhomogeneity, rather than modified gravity or dark energy. A perturbation is referred to as "subhorizon" if 
its wavelength is small compared to the Hubble length: X/Lh « 1. It has been suggested that small-scale density 
perturbations could cause the appearance of accelerated expansion without the need to introduce any form of dark 
energy, which is an appealing prospect 0, [H, 0, [H • The fact that inhomogeneity can systematically modify our inter- 
pretation of cosmological measurements was first realized by Ellis, who called it the "fitting problem" [H,[l3]. The 
basic idea is this: Due to the nonlinearity of the Einstein equation, the operators for taking spatial averages and for 
time evolution do not commute. This means that, although our universe is homogeneous in the mean, it will likely 
not have the same time evolution as that of the corresponding homogeneous universe. Nevertheless, we routinely fit 
distance data to FRW models, a procedure that introduces errors into the inferred properties of our Universe, and 
these errors will be present even for very large samples of Type la supernovae. 

Our goal in this paper is to calculate the lowest order fitting effect by calculating the cosmological constant density 
Da that one would deduce from a perturbed luminosity distance-redshift relation Dr,(z). If we treat cosmological 
fluctuations perturbatively and as a random process as suggested by the "fair sample hypothesis" [ll[ , then this fitting 
effect should be fundamentally nonlinear in the density contrast 5 = (p — (p))/ (p) , requiring that we work to at least 
second order in S. This is because the ensemble averages of first order quantities vanish. We model observations out 
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FIG. 1: The best-fit cosmological constant density Ha plotted as a function of the maximum redshift z ma x, for the choices 
Zmin = 0.01, 0.02, and 0.03. The horizontal dash-dot line shows the actual model value J7a = 0. 



to some moderate redshift z max ~ 0.1 -C 1. Within the corresponding comoving spherical region, the Hubble flow 
velocity vh is bounded above by 



< 



0.1 , 



(1.1) 



allowing us to use post-Newtonian expansions. There are two different velocity scales that occur, the Hubble flow 
velocity vh and the peculiar velocity v p . The corresponding dimcnsionlcss small parameters are 
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where A c ~ 10 Mpc is the wavelength of the dominant perturbation mode. In our computation, we will treat both of 
these parameters as being of formally the same order, and we will denote both by "e" for book keeping purposes. At the 
end of our computation we can identify terms that scale as s^e™ for different values of m and n. As mentioned above, 
we also expand separately in the fractional density perturbation 5. We will compute redshifts z(X) and luminosity 
distances HqD^(X) as functions of the affinc parameter A to third order in e and to second order in S. Combining 
these results to eliminate A will yield Dj, as a function of z. 

Using this expansion method, we find that the lowest order inhomogeneity-induced correction to the luminosity 
distance scales as \ADl\/Dl ~ S 2 (Hq\ c /c) ~ 10~ 5 . We then fit this relation to what one would expect from a 
homogeneous cosmological model which contains dust with a density Qm and a cosmological constant with a density 
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by maximizing a likelihood function. Here F(u) = u for a flat universe, T(u) — sinh(u) for an open universe, and 
F{u) = sin(u) for a closed universe. We find that the result for the cosmological constant density is dependent on 
the size of the redshift range for which we have supernova data. These results are summarized in Figure [T] For data 



from z r , 



0.02 out to a limiting redshift 



0.15, we find that the best-fit cosmological constant density is 
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57a ~ 0.004, and fl\ tends to get larger as z max gets smaller. The best-fit |Oa| also becomes larger as z m i n becomes 
smaller, since \AD^\/ becomes large on small scales. Although this ensemble averaged result is still quite small, 
we find that the variance can be a\ ~ 1 for a sample of 100 supernovae out to a redshift z max ~ 0.2. One implication 
of these results is that precision measurements of the cosmological constant from nearby supernova data require that 
we measure Dl(z) over a large enough redshift range, with a large enough sample. One could also try to correct 
for some of the effects of inhomogeneity, using available information about large scale structure and about our own 
peculiar velocity [H, Ell ■ 

The analysis presented here is more realistic than similar analyses within the context of simplified models of 
structure formation, such as the spherically symmetric Lemaitre-Tolman-Bondi (LTB) models [H], [H, [TtJ , Swiss 
cheese models [l8[ and their variants [H, 0, HH H2] . This is because we look at the full three dimensional problem, and 
assume that there are no bulk flows on cosmological length scales. There have also been analyses of the perturbations 
to the luminosity distance-redshift relation that go to Newtonian order [23|, [U HH, [26[ , that only consider super horizon 
perturbation modes [27],[28|], and that use Taylor expansions of the luminosity distance (2j|, which are most appropriate 
for long-wavelength perturbations. In contrast, we go to post-Newtonian order, we only consider subhorizon modes, 
and we fit to FRW models, so that we may fully address the "fitting problem" . 

Our analysis is also fundamentally different from those in Refs. [a, @, 0, H, HH , as we choose a different method 
for obtaining averaged expansion parameters. These authors average the expansion rate over a constant time slice, 
whereas we choose to calculate only observable quantities, namely the luminosity distance and the redshift, along the 
past light cone of the observer. We then combine these expressions into D^(z, 6, <fi), average over viewing angles and 
ensemble average, and then fit the results to what one would expect in a homogeneous model containing dust and a 
cosmological constant to find the best-fit value for fi^- This approach better simulates the process of gathering and 
analyzing supernova data, and it leads to a different result with a stable perturbative expansion. 

Refs. [a, [a0;ll] base their characterization of the expansion rate of the Universe on quantities that are not related 
to how observers have deduced the existence of dark energy. In these papers, perturbations are spatially averaged 
over a constant time slice. Such a spatial average is somewhat arbitrary, as it is dependent on the choice of spatial 
hypersurface. This is in contrast to the observable significance of D^(z). Refs. 0, H[ also use the synchronous 
gauge for their calculations, wherein there are metric perturbations of order S. Since S > 1 on small scales, this gauge 
is particularly ill suited to perturbation theory. In contrast, in Newtonian-type gauges the metric perturbation is of 
order S(Hq\ c /c) 2 <C 6. We explore this difference in Section VI. 

The organization of the paper is as follows: In Section II below, we introduce our coordinate choice, wherein 
we recast the Friedmann-Robertson- Walker (FRW) metric as an expansion around flat space, and in Section III we 
present the fundamental post-Newtonian optics equations that we will need for this calculation. We then explain our 
method of computation and calculate the necessary unperturbed quantities in Section IV. Here we also compute the 
luminosity distances and redshifts for a perturbed matter dominated universe, finding z and HqDl to second order in 
8 and to third order in e, and we find that we may write the lowest order correction to Dl(z) in terms of the peculiar 
velocity field. We then fit to a homogeneous model in Section V to find the best-fit J7a and its variance. The detailed 
redshift and luminosity distance equations are in Appendix A, the necessary results of second order perturbation 
theory are reviewed in Appendix B, and the averaging is discussed in Appendix C. Then, in Section VI we discuss the 
previous results in the synchronous gauge and show that one can choose coordinates and a definition of "acceleration" 
such that it appears as though there could be a larger fitting effect. We argue that such a result would be unphysical. 
A detailed discussion of transforming to synchronous coordinates is given in Appendix D. Finally, in Section VII we 
make our concluding remarks. As usual, Greek indices will be summed over all four spacetime dimensions while Latin 
indices will be summed only over the three spatial dimensions. We will also write 3- vectors in boldface and put arrows 
over 4- vectors. 



In general, certain coordinate choices allow us to conveniently recast the local metric as an expansion around flat 
space, as was first emphasized for the FRW metric by Peebles [33]. We will take advantage of such an expansion so 
that we may use the standard post-Newtonian formalism for this calculation. Starting with the usual FRW metric 
with c = G = 1 , 



II. POST-NEWTONIAN EXPANSION OF THE LOCAL FRW METRIC 
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we can define the new radial coordinate 



r = o(r)x 
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so that the line element becomes 



ds 2 = - (I - H 2 f 2 \ dr 2 - 2Hfdrdf + df 2 + f 2 dVL 2 



(2.3) 



where the Hubble parameter of a flat and dust-dominated FRW universe is H(t) — (l/a)(da/dr) = 2/3t; we will 
specialize to this case for the remainder of this paper. Now we change coordinates to the standard post-Newtonian 
gauge. In this gauge, the metric to first post-Newtonian order can be written as 



ds 2 = g^dx^dx" = - (l + 2$ + 2$ 2 ) dt 2 + 2Qdx l dt + (1 - 2$) "f l3 dx l dx 3 , 



(2.4) 



where jij is a flat spatial metric, the potential $ contains both Newtonian and post-Newtonian pieces, Q is the usual 
gravitomagnetic potential, and 

3$ + V-C = (2.5) 
is the gauge condition. Achieving this form for the metric entails transforming from t and f to t and r, defined by 



t 



1 



r 

3*2 



30i 4 



O 



(2.6) 



and 



t 2 „ ( r 4 
9P + \¥ 



(2.7) 



Then the line element becomes 
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to the necessary order in r. This metric is of the post-Newtonian form (|2.4[) if we define 
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Here subscripts "(0)" denote unperturbed, background quantities; we will add cosmological perturbations in subse- 
quent sections. The unperturbed density in the new coordinates is 
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and the continuity equation tells us that the unperturbed 3- velocity must be of the form vyq) = VfQ\d/dr, where 
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and where U(q) — dr/dt. Thus, we see that counting orders of e ~ v/c is equivalent to counting orders of r/t in these 
coordinates. Our coordinate choice and expansion method also have the consequence that the analysis of this paper 
is only valid for small redshifts. 

In general in the standard post-Newtonian gauge, the connection coefficients are 



r I = * , 



(2-13) 



(2-14) 



rf 3 - = -*ru - Cm 
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(2.16) 



and 
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to the necessary order in e, where T l - k is the connection associated with the flat spatial metric jij, which we will choose 
to be that of standard spherical coordinates (r, 9, </>), as in Ref. [33(. Vertical bars represent covariant derivatives with 
respect to 7^ . We will also need the Ricci tensor components 



and 



Rtt = V 2 $ 



Rij = V 2 $% 



Furthermore, the first post-Newtonian hydrodynamic and Einstein equations are 
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and 



~dt 



(v • V) v = -V ($ + 2$ 2 ) - C-(VxC)xv + 3$v + 4v (v • V) $ - u 2 V$ 



V 2 $ = Airp (1 + 2v 2 - 2$) , 



V 2 C = 167T/7V + V$ , 

in this gauge. The 3-velocity v is related to the 4-velocity u of the fluid by 

u= (u^u 1 ) =j(l,v l ) , 

where demanding that u ■ u = — 1 yields 

7 2 = 1 + v 2 - 2$ + 2$ 2 - 6$u 2 + v 4 + 2( ■ v . 
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III. COMPUTATION OF LUMINOSITY DISTANCE AND REDSHIFT 



A. Computing Dl(z) in a general spacetime 

In this section we will review how to compute luminosity distances and redshifts in a general spacetime, as in Refs. 
[33l [3^ |. Our analysis is initially similar to that of Ref. (29|, although they eventually rely on Taylor expansions 
around the observer's location. Such expansions are sensible for long-wavelength perturbations, but not for the short- 
wavelength perturbations that we consider here. We focus attention on a particular observer at some event V . In 
our application to perturbed FRW spacetimes, this observer will be at r = and at t = to for some fixed to. We 
consider the congruence of geodesies forming this observer's past light cone. Given the connection, we then find ray 
trajectories x a (X) by noting that the 4-momentum is k a = dx a /d\, and by using the geodesic equation 



(3.1) 



where we have defined d/dX = k a d a . Here the affine parameter A is chosen such that A = at the observer and 
A = A s < at the source. We also note that the 4-momentum is null. 
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The expansion 9 of the congruence of null rays is related to the area -4(A) of a bundle of rays by 



A dX 



We can find 9 by using the Raychadhuri equation 
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where we have defined the shear of the congruence 
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and where we require 9 ~ 2/ A as A 
given by the differential equation 



0, so that the area of the beam goes to zero at A = 0. The shear a 
da 



dX 



= -a9 + C a i3^k a k v i^ 



(3.5) 



where Cap^ is the Weyl tensor, and we have defined a null Newman-Penrose tetrad composed of the real 4-vectors 
k^ and m/*, and the complex conjugate 4-vectors £ M and t M . These satisfy the orthogonality conditions 



fc M m M = t% = 1 



and 



k^k^ — m^m^ 



ft. 



k — m^t^ 



, 



(3.6) 
(3.7) 



as in [29]. They are chosen at the observer and then extended along each geodesic in the congruence by parallel 
transport. We also choose the initial condition a = at A = 0. 

Once we find 9, we then find the luminosity distance as a function of the affine parameter at the source, 



D L (X S ) = lim 
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-X s (l + z) exp 
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(3.8) 



where AA corresponds to the size of the observer's telescope, which we set to zero. The right hand side of Eq. 
has a well defined, finite, limit as AA — ► due to the aforementioned initial condition placed on 9. Note also that the 
right hand side has an overall minus sign due to our convention that the affine parameter is negative. 
The rcdshift observed at A = 0, of the light emitted from the source at A s , is 
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where 
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u a k a = 7 {guk 1 + g ti k l + guv'k* + g io v r k j ) 



(3.9) 



(3.10) 



and where the subscript "s" will in general denote quantities evaluated at the source at the emission time and the 
subscript "o" will denote quantities evaluated at the observer at the observation time. By combining Eqs. (|3.8[) and 
(|3.9p we can, in principle, compute Dl as a function of z in a general spacetime. 



B. Computing Dl(z) to first post-Newtonian order 



Now we specialize the results of the preceding subsection to a perturbed FPW metric in the post-Newtonian gauge 
(12. 4|) . Our goal is to find both HqDl and z to order e 3 . At the observer, we have chosen r = and t = to and we have 
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normalized the 4-momentum such that k r = — 1. This implies that A f=a — r and r/t ~ — A/f ~ e to lowest order. We 
will thus need to find the right hand side of Eq. (|3.8|) to order Ae 2 so that we may find H Q D L to order e 3 . Because of 
this, we see that we will need the integral in the exponential to order e 2 , and therefore we will need to find X9 to order 
e 2 . Similarly, inspection of Eq. Q3.9P tells us to what post-Newtonian order we will need to compute the components 
of k a . To lowest order, g tt ~ 1, g ti = gn ~ e 3 , gij ~ 1, 7 ~ 1, and v l ~ e, and therefore we will need fc* to order e 3 
and we will need the spatial components k l to order e 2 . 

The post-Newtonian pieces of k a must be as small or smaller than order e 2 , as can be seen by noting that XT" V ~ e 2 
in the null geodesic equation (|3 . 1 1) . Given this assumption and the normalization of k a , Eq. (|3.ip reduces to 

Plugging in the connection coefficients from Eq. (|2.13jl - (|2.18j) . we find 
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Using the specified initial conditions, the solutions to these equations are 



£(r**) = --^*,* + 0(V) • (3.15) 



jfe* = 1 - 2$ - Cr - 2 / $dr' + O (e 4 ) , (3.16) 







k r = -1 + (e 3 ) , (3.17) 



k e = - r^e + Oie 3 ) , (3.18) 



^ = — K- n r^$ + O(£ 3 ) ; (3.19) 



and 



r sin U Jo r 

the integrals above are performed along the unperturbed ray, where t(X) = to + X and r(A) = —A. We can then find 
the perturbed ray trajectory by integrating Eqs. (|3.16[) - (|3.19p with respect to A. Most notably, Eq. (|3. 1T|) leads to 
A = — r + 0(Xe 3 ). This means that we can easily rewrite Eq. (|3.8p in terms of the radial coordinate r of the source: 



Dl = t (1 + z) 2 exp 







Ha* 









0(re 3 ) . (3.20) 



In order to find the expansion 8, we first need to find the shear, given by Eq. (|3.5p . The solution to this equation is 
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{X')C a(}lMV k a k v tn f3 dX'; (3.21) 



since |fc| ~ |t| ~ 1, it turns out that the lowest order shear is a ~ e 2 /A. Inserting \a\ 2 ~ e 4 /A 2 into the Raychaudhuri 
equation (|3.3p gives a contribution of order e 4 /A to the expansion 9. However, we already know that we only need 9 



8 



to order e 2 /A, and so this contribution is negligible for our purposes here. Neglecting shear and defining 59 = 9 — 2/A, 
we rewrite Eq. (|3.3[) as 
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* / (r') 2 V 2 $dr' + 0(^-) , (3.23) 



The solution to this is 



where we are using A = -r + 0(Ae 3 ). Using this result in Eq. (t3T20|) yields our final result for the post-Newtonian 
luminosity distance 
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We now turn to evaluating the redshift z as a function of the affine parameter A. Equation (|3.9p is the general 
expression for the redshift, and it depends on u a k a at the source and at the observer. To order e 3 , using Eqs. (|2.4p . 
(I2.26p . and our solutions for k a , we obtain 

u a k a = g al 3U a kP 

i i r ■ 

= -l-v r - -v 2 + $ + 3v r <f>- -v r v 2 + 2 / <I>dr' + v e k e + + O (e A ) , (3.25) 

where k 9 and fc^ are given by Eqs. (|3.18p and (|3.19[) . respectively. Therefore, the post-Newtonian redshift is 

(u a k a ) s 
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= 1 + v r s - v r a + $ - $ s + \ (v 2 s - v 2 ) + (v r Q ) 2 - v r v r s - 2 f $dr' + (vgk e + v^k*) - (v e k e + v^) s 

2 Jo 

+ + + d> < - 3$ s < - ^ 2 « - <) + «) 2 K - <) + ^ « - <) + O (e 4 ) ■ (3-26) 

In Eqs. (|3.24p and (|3.26p . the right hand sides are evaluated at r = —A and t = tg + A. Recall that subscripts "o" 
denote quantities evaluated at the observer where r — and t — to, while subscripts "s" denote quantities evaluated 
at the source (t(X), r(X), 9, (j)). 

IV. ADDING DENSITY PERTURBATIONS 
A. Basic method 

In this section we apply the formalism of Section III to a spherical region in a perturbed FRW spacetime. We will 
describe that region using the post-Newtonian metric (|2.4p . We expand the metric functions <I> and Q and the fluid 
3-velocity v* in powers of the density contrast 5 as 

$ = * (0) + + $ (2) + O (<5 3 ) , (4.1) 

Ci = Ci(0)+Ci(i)+Ci(2)+O(S 3 ) , (4.2) 

and 

^=4 ) )+4)+4 !) + 0(£ 3 ) , (4.3) 
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respectively. We also expand the null geodesic x a and 4-momentum k a = dx a /d\ as 
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respectively. For the remainder of the paper, quantities that are zeroth order in S will be denoted by a subscript "(0)" , 
first order by a subscript "(1)", and second order by a subscript "(2)". Also henceforth "first order" and "second 
order" will always refer to orders in 6, not e, unless otherwise specified. 

In the perturbed spacetime, we will calculate the redshift z and luminosity distance Dl as functions of the obser- 
vation time to, of the affine parameter A along the past-directed null geodesic, and of the 4-momentum k of photons 
at r — and t — to- We parameterize this future-directed null vector k in terms of angles 9 and cj), in such a way that 
k r = — 1 and k is in the direction (9, <ft) at r = 0. We can thus express Dl and z as functions of A, 9, and 4> at fixed 
to, and by eliminating the affine parameter A we can compute Dl(z, 9, 0). 

We can then take an average over angles to find Dl(z), where we must take some care since there are two sets of 
relevant angles. There are the angles {9, <f>) which parameterize the direction of k in the observer's rest frame, and 
then there are the coordinate angles (9,(f>). We will need to average over (9,<fi). This means that we will need to 
know the relationship between the related infinitesimal solid angles dil 2 and dfl 2 . We define Cartesian coordinates 
(x ,x ,x ) in terms of the polar coordinates (r,9,(f>) in the standard way. An orthonormal set of basis vectors for 
the observer's local Lorentz frame can be obtained by renormalizing the coordinate basis vectors d/dt and d/dx l and 
boosting. The result is 
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dxi 



The angles (0,4>) are defined by 



with 



k = k 8t 



dx 1 



n = (sin 9 cos (j>, sin 9 sin 4>, cos ( 
while the observer's angles (9, </>) are defined by 

k cx &t — fi €i , 

with 



n = ^si 



sin a cos 0, sin o sin 0, cos ( 



By inserting (|4.6[l and (|4.7|) into (|4.10|) and then comparing with (|4.8jl , we find 



(4.6) 



(4.7) 



(4.8) 
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This gives 



n cx n + $ c n - v + - (v • n) v + O (e 3 ) 



d 2 n = d 2 n [1 - 2 (v • n) + O (e 2 



(4.12) 



(4.13) 



After averaging over viewing angles, we find the expected value of Dl(z) by taking an ensemble average, wherein 
we treat the density perturbation 6 at any fixed time as a homogeneous random process. Once we have the averaged 
Dl{z), we can then analyze these data in terms of a homogeneous universe to see if we would find an apparent 
acceleration. Assuming Gaussian uncertainties, we perform a chi-squared fit to a FRW model with a matter density 
CIm and a cosmological constant density Qa- 
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B. Unperturbed quantities 

In the unperturbed background, everything is spherically symmetric, and the line element in our coordinates is 
given by Eq. (|2.8p . The background four-momentum k$ is purely in the t — r plane, and is given by Eqs. (|3.16p and 
(13T7) to be 



fc( 0) M) = i 



9t 2 135i 3 



O 



t 4 



(4.14) 



and 



fcfo(r,t) = -l- 



4r 3 / r 4 

27*3 VF 



(4.15) 



Since k*, s = dt/dX and fc^ = dr/dX, we can integrate and invert these equations to find the unperturbed ray 
trajectory; keeping in mind the conditions that r = X = and t = to at the observer, we find 



and 



t(X) = t + X 



r(A) = -A 



2A 2 






-(f)] 







1 + 



(4.16) 



(4.17) 



in the unperturbed background. 

Using this, we can use the solution (|3.23p to the Raychaudhuri equation to find the background expansion 0(o), 



'(O) 



(A) 



A 9t 2 



X + O 



A 



(4.18) 



Then the zeroth-order luminosity distance is given by Eq. (|3.24|) to be 



D 



L(Q) 



(1 + z) r I 1- 



r 8r 2 



1 -i + w +0 \^ 



where we have defined Ho = 2/3io- The zeroth-order redshift is found from Eq. (]3 . 2G[) , 

2r r 2 4r 3 /r 4 ' 

and we eventually find the expected Dl(q\{z) by inverting Eq. (|4.20[) and plugging the result into Eq. (|4. 19[) : 



(4.19) 



(4.20) 



D m {z) = jr 

no 



l + \ z -\ z 2 + 0^) 



(4.21) 



Thus, for the background, the best-fit cosmological constant density is J7a = and the deceleration parameter is 
qo = 1/2. 



C. Second order perturbed optics 

The perturbed post-l-Newtonian line element is, from Eq. 



ds 2 



1 + 2$ 



(0) 



2$ { i) + 4$ { o)$ { i) + 2$ (2) + 2$f 1} + 4$ (0 )$(2) dt 



+ 2 (Ci( ) + Ci(i) + 0(2)) daJ*tft + (1 - 2* (0) - 2* (1) - 2* (2) ) i^dxUx 1 



(4.22) 
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and the perturbed luminosity distance (|3.24[) is denned to be 

D L = (1 + z) 2 E L , 

where 

E L {r,6,(ft) = E Lm {r,9^)+E L(1) {r,9,cft)+E L{2) {r,9^)+0{8 :i ) 



(4.23) 



* (2) ) dr" 



O (<5 3 ) 



(4.24) 



Jo r Jo 

and where we have pulled out the factor of (1 + z) 2 for simplicity. We then find that the order S perturbation is 



and the order 5 2 perturbation is 



Eh(2) — 



2 r 
3H t 



Jo > Jo 



(4.26) 



where Hq = 2/3io- In general, all of the terms involving potentials and velocities in these equations, and in those 
that follow, are evaluated along the zeroth-order, unperturbed, geodesic. 
We can now calculate the perturbed redshift 



z(r, 9, (ft) = z (0) (r, 9, (ft) + z (1) (r, 9, (ft) + z (2) (r, 9, (ft) + 0(6 3 ) 
from Eq. (|3.26|) . using our knowledge of the zeroth-order quantities, to find 



(4.27) 



Z{1) 



•Mi) 



Mi) 



$o(i) - *«( 



1) + § (<(1) " <(!)) - 2 l ®(¥ r ' + I , 



and 



: (2) 



V s(2) ~ V o(2) 



r r n 



2r r 



3f 2 



(4.28) 



2r 
3t 



M2) + — (<(2)-<(2) 



2r 



^(^(1) -<!)) + (<(1)) 
(<(!)) +(<(! 



2r 
3t 



^(l) U o(l) 



where the first order perturbation to the null geodesic is 



x\i) = ~ [ Hi)dr' 
Jo 



(4.29) 



(4.30) 



All of the quantities above are evaluated along the zeroth-order geodesic, and the integrals are performed along an 
unperturbed central ray where r(A) = —A and i(A) = t + A. 

Now we have found the redshift z and luminosity distance HqD^ as functions of affine parameter A and initial 
4- momentum k a , to second order in S and to third order in e. Adding the redshift equations (|4.20p . (|4.28|) , and (|4.29D 
yields z(X,9,<ft). Similarly, the luminosity distance Dl(X,6,^>) is found from adding Eqs. (|4.21|) . (|4.25[) . and (14.26[) . 
after replacing the factors of (1 + z) 2 . Inverting z(X, 6, (ft) perturbatively, in terms of either S or e, gives us A as a 
function of z. Plugging this into Al(A, 6, (ft) yields an expression for Dl(z,6,<P). We then angle average this and 
then ensemble average, assuming that density fluctuations at a given cosmic time are a homogeneous random process. 
Details of this full procedure are given in Appendices A, B, and C, and the result is 



D L (z) = — 1 + -z- -z 2 
V ' H a \ 4 8 



AD L (z) 



(4.31) 
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where ADl(z) depends on the two point correlation function. We will only need the lowest order piece of this, which 



is 

4 ^> = -3V(A)<"*>> +0 (t?) ' ( " 2) 

using z « Hqt to lowest order. The function / is related to the velocity two point correlation function (see Appendix 
C): 

t , \ 3(n- v(r ,t)n- v(r +rn, t)) 

m = ¥i^W) 1 ' (4 - 33) 

where n is a unit vector that defines the viewing direction and ro is an arbitrary location in space. Note that f(r) 
is independent of time, even though (v(ro,£) 2 ) does depend on time. This is because the time dependences of the 
numerator and denominator cancel. 

The perturbation to the luminosity distance is proportional to 



2 



4 



((V$ o( i)) 2 ) ; (4.34) 



this qualitative scaling has been argued for in Refs. *35| and .'Hi . We can Fourier transform $(i), in terms of a 
wavevector k l (not to be confused with the previously-defined 4- momentum) Q, 

so that we may write the average of (V<£>(i)) 2 as a sum over modes: 

<(V$(1))2H 4^ 4 /~F A2(fc) ' (4 ' 36) 

where A(k) is the dimensionless power spectrum of matter density fluctuations at the present time, defined by 

/oo 

d(\nk)A(k) 2 . (4.37) 
-OO 

We adopt the following power spectrum 

A s W =C*(A)V(A), (4 . 38) 

where the factor of (k/Ho) 4 reflects a Harrison-Zel'dovich flat spectrum, the amplitude C = 1.9 x 10~ 5 is set by 
observations, and T(y) is the transfer function. The BBKS transfer function |37j is a good fit for T in the absence of 
dark energy, 



T(y) = lni \\ 2 . My) \l + 3.89y + (16.1y) 2 + (5.46y) 3 + (6.71y) 4 j ^ , (4.39) 
2.34y L J 



where 



k ^ IA AM 

y = - — = r . (4.40) 

k eq f^h 2 MpC _1 

Here we show the most general form of the transfer function, where 9 = per/I-G&Pj (not to be confused with the 
expansion 6) is the density of relativistic particles divided by the density of photons, Qx is the density of cold dark 
matter, and h = i?o/(100 km s _1 Mpc -1 ). We choose fix = 1 for our analysis. 
Using this spectrum, 

(( V$ o(i)) 2 ) = — / ydyT 2 (y) , (4.41) 
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FIG. 2: The function 1 + f(r) plotted versus k eq r, where k eq is the wavenumber of the dominant perturbation mode. 



where k eq = 1/A C = ^ h^J^ Mpc" 1 » 3OOOfixh0- 
transfer function in Eq. (|4.39|) . So we finally find 



1 / 2 i?o and the integral is approximately 2.31 x 10 2 , using the 



((V4 



o(l)J 



9 x lO"^ 



0.27 



and therefore 



<t, 2 (1) ) « 3 x 10- 



A 



0.27 



h 

07 



h 

0/7 
n 2 



Using the power spectrum (|4.38p we also find 



(n • v(r , i)n • v(r + rn, £)) 



H 2 



ydyT 2 (y) 



3 J0 



3 J2 



(4.42) 



(4.43) 



(4.44) 



where jo and J2 are spherical Bessel functions of the first kind, defined in Eqs. (|C13|) and (|C14[) of Appendix C. We 
plot 1 + f(r), found by combining Eqs. (|4.33|) and (|4.44[) . in Figured Note that this becomes negative for k eq r > 10. 
Note also that we have not used any truncation of the power on scales that are nonlinear. If we instead were to impose 
a high-fc cutoff, so as not to include the effects of any modes that have A 2 (A;) > 1, then this would lead to differences 
of a factor of about two. A different approach would be to include the quasi-linear regime, with the power spectrum 
given from N-body simulations [38j |. 

We will specialize to k eq /Ho = 1000 for the rest of this paper, which yields 



8.34 x 10" 



(4.45) 



In Figure [31 we show how the perturbation ADl(z) scales relative to the unperturbed luminosity distance -Dl(o) (z), 
for the choice k eq /Ho = 10 3 . Note that we are plotting the logarithm of the absolute value, as the perturbation 
changes sign from positive to negative as one looks at larger distances. By inspection, it becomes clear that ADl(z) 
is not actually a perturbation for very small redshifts, i.e. for where \ ADl\/ D L ( Q j ~ 1, and thus our computation of 
ADl is no longer valid in that regime. Indeed, it is well known that the peculiar velocities of objects within the Local 
Supercluster are not small when compared to their redshifts. However, this will not be a problem in practice, as Type 
la supernovae at such small redshifts are typically not used for cosmological parameter fitting. We will eventually 
take this breakdown of perturbation theory into account by imposing a lower cutoff z m i n when we fit our data to a 
theoretical model. By eye, we see that it should be safe to choose z m i n ~ 0.01. 
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FIG. 3: The relative size of the perturbation \og\\ADh{z)\/ Dl(q)(z)] plotted versus log(z), assuming that the dominant per- 
turbation wavelength is 10 3 times smaller than the Hubble scale: k eq /Ho = 10 3 . 



V. THE PERTURBATION TO THE INFERRED COSMOLOGICAL CONSTANT 

A. Finding the best-fit FRW model 

We may now find the inferred cosmological constant and deceleration parameter by analyzing Eq. (|4.3ip within the 
context of what one would expect in a homogeneous model. The lowest order perturbation to the luminosity distance 
depends on the difference between the peculiar velocities at the source and at the observer, and so the question that 
we now ask is: How do peculiar velocities and their correlations affect inferences drawn from data about cosmological 
models? We cannot simply Taylor expand Eq. (14.31|) around the observer to find qo- This is because / varies on 
short lengthscales of order k~ q ^10 Mpc, so that a Taylor series expansion would effectively mean computing go from 
Dl(z) within this unrealistically short lengthscale. A good alternative then is to fit the perturbed luminosity distance 
over a finite range of rcdshifts to what one would expect in a homogeneous model with matter and a cosmological 
constant. 

Suppose that the observer can measure redshifts {zi} for a set of distant objects arbitrarily well. From the distance 
determinations {D^}, the observer can compute {r^ = D^i/^l + ^i)}, and we can therefore take {zi,ri} to be the 
data gathered by the observer. Suppose also that in actuality the Universe is spatially flat with Hubble parameter 
Hq and matter only. Let 

n = H^ 1 [F{z l ) + A l {z i )] (5.1) 
be the physical value of r^, where for a flat matter-only cosmology 

and Ai(zi) (not to be confused with the matter perturbation power spectrum) is the non-FRW contribution to r^, 
from fluctuations via velocity differences. From Eq. (|4.32p . we find the ensemble averaged perturbation 

ycoa(k eq z t y/H ) _ ^ si^k^y / 'Hp) _ 6 cos(fc eg z t y/iJ ) 
k eq Zi/H (k eq Zi/H ) 2 (k eq z t /H ) 3 y 




(5.2) 



(^21,3 f oo 

A,(z,) « — $L / dyT \y) 
n o Jo 



sin(fc eg z i y/iJ ) 
{kegZi/Ho^y 2 



(5.3) 
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The observer fits the data to a FRW model that is slightly curved and has a small cosmological constant. The 
fitted model is then 



dz 



o H(z) 6 \J H{z) 



dz 



(5.4) 



where k = (Am + — l)-ff? if and 



H 2 (z) = H% t [Q M (1 + zf + (1 - Q M - fi A )(l + z) 2 + 



H%{l + zf 



1 



(l-ft M )z ft A z(2 + z) 



1 + z 



(1 + *) 3 



(5.5) 



here Hfn is the fitted Hubble parameter, and Q,m and f2 A are the density parameters for matter and for the cosmo- 
logical constant, respectively. Let us work to first order in 1 — Qm and A , a simplification which ought to suffice as 
long as Aj <C 1. Thus, the fitted model is 



H lt = Hy l 1 t [F{z l ) + {i-n M )G{z l ) + n A i(z i )} 

= Hj* [F( Zi ) +e M G( Zl ) + e A I( Zi )} , 
where F(zi) is the same as before, and we have defined 



G(z) 



dz z 



2 J (l + z)5/2 6 



[F{z) 



and 



I(z 



dzz(2 + z) 1 



(1 + z) 9 / 2 6 



(5.6) 



(5.7) 



(5.8) 



There are three fitting parameters: Hfu, tu = 1 — and e A = A . 

From the data and our model we can compute a likelihood function. Assuming Gaussian uncertainties this will be 
the exponential of 



x 2 



[n - r{ lt {zi) 



-- E 



1 . ^ [(#0 1 _ H fit) F i + H o lA '' ~ H fit( € MGi + e A Ii) 
= ~2^ a} ; 



(5.9) 



where Gi is the estimated uncertainty in the value of ri inferred from observations and Qi = Q(zi) for Q = F,G, I. 

The next step is to maximize x 2 with respect to the parameters of the fit, which will lead to a set of coupled 
nonlinear equations. To simplify, let us linearize in the small parameters 6m, £a, {Ai} and h — Hf it /H — 1. The 
resulting equations are 



(Aj-Fj) - e M {G i F i ) + e A (7 i F i ) - h{F?) , 



(5.10) 



(Aid) = e M {Gl) + e A (IiGi) - h(GiFi) , 



and 



(AJi) = e M {I t G t ) + e A (lf) - h^F,) , 
where we have defined the average (Qi) = J2i Qi/(N<jf). Solving for the parameters of the fit, we get 



(5.11) 
(5.12) 



v- 1 



(AifiX^fiXG?) - (IiGJiGiFi)) + {A t G i ){(UG l ){Ff) - (/^(G^)) 



+ <A^)((G^> 2 -(G 2 }<F 2 )) 



(5.13) 
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e M = V- 1 {/\ l F l ){{G l F l ){lf)-{I l F l ){I l G l )) + {^G i WiFi) 2 -{lf){Ff)) 



-(A i 7 i )«7 i G i )(F?) - (GiFJihFi)) 



(5.14) 



and 



h = -V- 1 



(A^H^G,) 2 - (7 2 )(G 2 )) + (A i G i )«G i F i )(7?) - 



-(A I 7 i )((7 I F,)(G 4 2 ) - (GiFJilid)) 



(5.15) 



where 



V = (IiG t f(Ff) ~ 2{I i F i ){I i G i ){G i F i ) - (7 2 )(7; 2 )(G 2 ) + (7 2 )(G lJ F i ) 2 + (7 i F l ) 2 (G 2 ) 



(5.16) 



These are fairly general for small Aj, and show that there may be contributions to e\, 6m, and h from velocity 
fluctuations. 

Next, we need to compute the averages. To do this, we recall that F corresponds to comoving radial coordinate, 
modulo a factor of H^ 1 . To the order of approximation underlying our calculations, we can take the comoving source 
density to be uniform. Moreover, we do not need to worry about Malmquist bias, at least for Type la supernovae, 
which are very bright. Let us also assume that all of the {of} are the same, to keep the problem as simple as possible. 
Then of drops out of our expressions for cm , e\, and h, although it remains in their uncertainties. We suppose that 
our source catalog extends to some maximum value F max , with a corresponding maximum redshift z max . It is worth 
remembering that F < 2 is an absolute upper bound, and that for z < 1, F < 2 — \pl ss 0.6, so we will be dealing 
with relatively small values of F typically. Moreover, as we have already noted in Figure [3l our small Aj assumption 
breaks down below a minimum redshift z m i n < 0.01, but this is not a problem as no supernovae below this redshift 
have ever been used for cosmological model fitting [3, ■ So we will assume a lower cutoff for all of our sums of F m i n . 
Then, for example, 



F 3 

max 



F 3 



dF F 2 F 2 



(5.17) 



and Eqs. (|5.2|> . (|5.7|) . (|5.8[) . and (|5.16[) give the lowest order result, assuming that F^ ax > F^ in , 



5268480 max ' 

Keeping only lowest order terms in F max in the numerators of Eqs. (|5.13p . (|5.14[) , and (|5.15[) as well, we get 



(5.18) 



5268480 
16 



784 F 4 28AF 15 ' ^d()F 6 

i <^t:j_ rnax max ^ max 



3(a. ( 7; 2 



(5.19) 



€M ~ — 2eA , 



(5.20) 



and 



5268480 
16 



(AiF) (A t F 2 ) , SiAtF 3 ) 



4487 12 1687 13 784F 4 

^ max max ' w max 



(5.21) 



We see that if A, cx F i: then ej^ is zero, because the three terms in Eq. (I5.19[) cancel. This means that if Aj arises 
from velocity correlations, it is only the correlation function of velocities at two separated points that matters, not 
the RMS velocity at a point. Also note that, for this fitting procedure, the deceleration parameter is still go = 1/2, 



Aq Q = q - 
1 



1 



= 2^ 



(aa/H 2 ) 
1 



1 - 2ft A ) = 7t(-£a/ - 2e A ) = -(2e A - 2e A ) = 



(5.22) 
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from Eq. (|5~20j) . in agreement with Refs. [27j], 01, and 

The perturbation Aj, given in Eq. (|5.3[) . depends on the correlation function f(r), and so it does contribute to e\. 
For z m i n = 0.02 and z max = 0.15, we numerically integrate to find that the best-fit cosmological constant density is 
S!a ps 0.004. Table [T] gives a few more results for the best-fit values for eA, cm, and ft. as a function of the two limiting 
redshifts z m i n and z max in the continuum limit, where we have made the assumption that the number of sources N 
is very large: N — > oo. In this limit, Ai{z%) — > A(z) and 



e A = / dFw(F)A(F) , (5.23) 

where we have the weighting function 

= 5268480 / 3F 3F 2 | F 3 \ 

WM " 16 1 784^ 280F4 _ + 140F6 a J ' ( ' > 

We also plot these results in Figure [TJ in the Introduction. Note that 17 a may be positive or negative, depending on 
the redshift range, since AD^ changes sign in the region of interest. 



TABLE I: Best-fit parameters in the continuum limit for a few values of the source catalog limiting redshifts z m in and z, 
also for the choice that the dominant perturbation wavelength is 10 3 times smaller than the Hubble scale: k eq /Ho = 10 3 . 



Zmin 


Z-m a x 




I-Om 


H fu /Ho - 1 


0.01 


0.1 


-0.018 


0.036 


-4.3 x 10" 5 




0.2 


0.0016 


-0.0032 


4.0 x 10 -5 


0.03 


0.1 


0.0037 


-0.0074 


7.1 x 10~ b 




0.2 


0.0020 


-0.0040 


4.7 x 10 -5 



In order to test the robustness of these continuum limit calculations, we have also applied our fitting procedure to 
randomly-generated catalogs of synthetic redshift data. To generate a data point Fi for such a catalog, we assume 
that the quantity {Ff — F^ in )/ {F^ ax — F^ in ) is distributed uniformly between and 1. In this way, we create catalogs 
of N = 100 data points, wherein each data point is a value of Fi for a source with a random location. For each data 
point, we use the ensemble averaged formula for ADl(z) to find Aj. We then fit these data to a homogeneous model 
as outlined above, using sums instead of integrals. Using 20 randomly-generated catalogs, the average best-fit values 
for fi A are summarized in Table [Til along with their standard deviations. We also found the best-fit cosmological 
constant with 50 catalogs for z m i n = 0.02 and z max — 0.15, to find Ha — 0.005 ± 0.001. 



TABLE II: Best-fit parameters for 20 catalogs of N=100 samples each, for a few values of the source catalog limiting redshifts 
Zmin and Zmax- We have also made the choice that the dominant perturbation wavelength is 10 3 times smaller than the Hubble 
scale: k eq /Ho = 10 3 . 



Zmin 


Zmax 




0.01 


0.1 


-0.020 ±0.002 




0.2 


0.002 ±0.001 


0.03 


0.1 


0.014 ± 0.001 




0.2 


0.0025 ± 0.0004 



B. Variance 



Although the best-fit values for Da of the previous subsection are very small, we must keep in mind that they 
are derived from the ensemble averaged perturbation to the luminosity distance. For a given source, this ensemble 
averaged perturbation will be far smaller than the leading order perturbation, which depends linearly on the peculiar 
velocity. This linear perturbation will be the main source of the variance in the best-fit parameters, and this variance 
should overwhelm the systematic error for typical supernova sample sizes. This complication was pointed out by Ref. 
[40| and it was shown to cause errors of A^a ~ —0.04 for a sample of actual nearby supernovae in Ref. fl2l |. 

Consider our expression for the best-fit in terms of N discrete sources, rewritten as a weighted sum, 



i 



(5.25) 
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What we have computed is the ensemble average of this, 

(°A> = ^$>(^)<Ai> . (5.26) 



N 

The variance is then 



4 



((n A -(n A )f) = (tf A ) + o(5 3 ) 

= ^E w WWfi){AA')- (5-27) 

which has two types of terms contributing: those with i = j and those with i ^ j. Separating these, we have 
u\ = erf + cr| , where 

^^E^XA, 2 ) (5.28) 

i 

and 

In the continuum limit N — > oo, the first piece of the variance becomes 

^ " ^ pi - /^^ F 2 dF W 2 (F){A 2 (F)) (5.30) 

max JQ 

where, from Eq. ()5.3|) . 

<A 2 (F)) = <A 2 (iJ r)) = (n • [v(r) - v(0)] n • [v(r) - v(0)]) ~ {v 2 D ) . (5.31) 

The integrand in Eq. (|5.30[) is integrable as F —> 0, and so the quantity o\ is to a good approximation independent 
of z m in for small z m i n . Thus we can for simplicity take z m i n — 0. After integrating, we find 

a 100 / (v 2 ) \ /w\~ 6 r _ 9 , 

ai ~ at (.^To^ JUtJ • (5 - 32) 

For a source catalog of 100 sources out to a limiting redshift z max = 0.2, we find that this variance is significant: 
a 2 ~l. 

The second piece (|5.29[) of the variance does not depend on the sample size, although it does depend on F max . In 
the continuum limit, 

f"" F 2 dFw(F) [ {F'f dF'w{F')(A(F)A{F')) (5.33) 
Jo Jo 



* P6 

max J 

where 



(A(F)A(F')) = \{v 2 o(l) ) 



f[---)-f(-)-f(- 



(5.34) 



Plugging Eq. (|5.34jl into Eq. (|5.33p . then using Eqs. (|4.33|) and (|4.44p . and then finally doing some rearranging, we 
find 



2 / 2469606*^0 
a 2 ~ 



\ F 5 k 



Jo y 



i 2,k e q Frnax 

l[ —Ho— y 



2 



(5.35) 



where 

I(q) = dx (^^x -^x 2 + ^x 3 ^j (sin qx-qx cos qx) . (5.36) 

This result for erf does not depend on the sample size, as it only depends on the size of the redshift range F max , 
making it a measure of cosmic variance. By integrating numerically, we find that it scales roughly as F~f [X and 



a 2 ~ 0.03 



(it)" ■ <" 7 » 



For comparison, Ref. [12j uses a sample of 115 supernovae up to a redshift z max — 1.01, and they find an error from 
the data of A£7a = —0.04. For this same scenario, we estimate |A57a| ~ 0.01, from the sum of Eqs. (|5.32p and (|5.37p . 
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VI. CONSISTENCY WITH PRIOR RESULTS 

The method of analysis that we have presented in the previous sections differs from that of Refs. [1, H, 0, H] . This 
is because of (i) a difference in gauge choice and (ii) a fundamental difference in the definition of what constitutes 
"acceleration" . We have chosen to use the standard post-Newtonian gauge, and to define acceleration as being based 
on fitting the luminosity distance-redshift relation to that of a homogeneous model containing dust and a cosmological 
constant. As this definition of acceleration is based only on observable quantities, performing our calculation in other 
gauges gives us the same results. 

In contrast, Refs. [H, @, 0, [1] calculate the cosmological expansion rate, averaged over a constant time slice. The 
motivation for doing this comes from the spatially-averaged Friedmann equations, also called the Buchert equations 
[4lj |. In particular, Ref. 0] defines the effective coarse-grained scale factor bd in terms of the average matter density: 
(p)d <x a D 3 , where the angle brackets ()r>, with subscript D, denote an average over a spatial hypersurface D at a 
given time. Then Ref. Q defines the coarse-grained Hubble rate 

Hd = — = \{6) D (6.1) 



and the effective deceleration parameter 



Hd 



1 . (6.2) 



These measures of acceleration are somewhat arbitrary since the deceleration parameter (|6.2p depends on the spatial 
hypersurface over which one averages. Refs. @, @, 0, Q use constant time slices in the comoving synchronous gauge. 
In this gauge, the perturbation to the the expansion 9 is related quite simply to the perturbations to the trace of the 
connection; from Ref. Q, 

{0 { i)) D = - a {T\ l{l) ) D , (6.3) 

and similarly for #/ 2 )- Ref- [61] claims that spatially averaged perturbations could become quite large, which implies 
that our perception of the expansion rate of the Universe is significantly affected by inhomogeneity. The culprit is 
the appearance of terms in T l ti ^ with large numbers of spatial gradients, which naturally appear in the synchronous 
gauge. These higher derivative terms, which do not appear in our method above, lead to a perturbative instability, 
wherein terms higher order in perturbation theory do not get smaller as expected. 

Although the results of the previous sections appear to differ from the claims of Refs. [!, H, 0, [1], m fact the large 
fitting effect claimed in those papers arises at a higher post-Newtonian order than we have computed. In this section 
we show that our results are consistent with theirs to the order we have computed. Our method of computation could 
be extended to higher post-Newtonian order, which would allow for a detailed confrontation with their claims. 

However, we believe that our result of a small fitting effect is robust, in the sense that it will not be altered by the 
inclusion of effects that are higher order in e and/or 5. This belief is based on the structure of the post-Newtonian 
expansion of Einstein's equations, and on the fact that we are computing a gauge-invariant observable. If this is true, 
then our conclusion is in disagreement with Refs. [H, SSSl- 

We believe the most likely reason for the disagreement is that we compute a gauge-invariant observable that is 
directly and uniquely related to supernova observations, whereas the quantities computed in Refs. 0, 0, 0, H| have 
some arbitrariness and are not directly related to observations. The proposal of Refs. 0, 0, 0, [f| that there might be 
a large backreaction effect in terms of qo does not necessarily imply that observers will measure large deviations from 
FRW dynamics. As mentioned above, spatially averaged perturbations are dependent on one's coordinate choice, 
in the sense that a constant time hypersurface in one coordinate system is most likely not going to be a constant 
time hypersurface in a different coordinate system. These averages are unlikely to be directly observable, and are 
not uniquely related to the cosmic acceleration inferred from cosmological observations. As Hirata and Seljak [28| 
remarked, we "cannot cover the entire universe with astronomers so as to measure spatially averaged quantities" such 
as Hp. It is possible that the measure of acceleration (|6.2p could be large while the observed acceleration is small. 

We now turn to showing consistency of our results with those of Refs. [a, H 0, [H to the order we have computed. We 
take our metric (|2.4p and transform it from the post-Newtonian gauge to the synchronous gauge. We then compute 
from the transformed metric the perturbation to the Hubble rate. The relative size of the difference between Hjj and 
the expected FRW value H determines whether or not there will be a large fitting effect. As an example, we will now 
compute the ratio 

H D -H _ AH _ (0(p +0( 2 ))z? 

H ~ H 3H [ ' ' 
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where the spatial average involves integrating with respect to the perturbed volume element dV = y/g spaced 3 x, where 
9spa.ce is the determinant of the spatial part of the metric. Note that the quantity that we define as AH/ H differs 
from what is computed in Refs. @, 0, @], although we do find the same qualitative result at the end of the 
day. Below wc show that this quantity is small to Newtonian order, in correspondence with what was found in [7|, 
even though it involves a sum of terms that can be large individually. The reason these terms are large is that in 
synchronous coordinates metric perturbations can be of order 5, which may be of considerable size even though there 
are no large gravitational potentials anywhere in the Universe. By contrast, in our calculation based on standard 
post-Newtonian coordinates, metric perturbations are at most of order e 2 <5, which is always small. In this sense, 
perturbation expansions are much better behaved in the standard post-Newtonian coordinates than in synchronous 
coordinates. 

We start by reviewing the transformation from standard post-Newtonian coordinates (|2.4p to synchronous coordi- 
nates; a detailed discussion is presented in Appendix D. Begin with the second order perturbed FRW metric in the 
gauge 



ds 2 = a 2 { V ) [- (1 + 2$ (1) + 2$ (2) ) dr, 2 + (1 - 2$ (1) - 2$ (2) ) S l3 dX l dX^] , 



(6.5) 



where we are now using conformal and Cartesian coordinates for simplicity, and we will only need to work to Newtonian 
order. We can then define the new coordinates r and x l by 

" (v«> (1) ) 2 



1 



3(1) 



1 



(2) 



45 



+ O (r e 4 ) + O(t 6 3 ) 



(6.6) 



and 



12o $ (D,<i*(i)j + ( 5 * e 



0(x l S 3 ' 



(6.7) 



where these potentials are fixed physical quantities, evaluated at (r, x l ), and these spatial derivatives are in terms of 
the new coordinates. We are also assuming that we have the growing mode only, for which we have the power law 
scalings $(1) oc r° and $(2) oc r 2 . Then the line element becomes, to lowest order in e, 

ds 2 = a 2 (r) [— dr 2 + g^dx 1 dx?~\ 



a 2 (r) -dr 2 + 



r 2 t 2 t 2t 4 



O e 2 



dx l dx j \ (,6.8) 



which is now in a synchronous gauge. Note that the metric now has perturbations of order e°S ~ S. These order 5 
perturbations will lead to the appearance of large terms in AH/H, which will cancel when averaged. Then we find 



a 3 (r) 



1- y V 2 $ (1) 



yj 9 'space 

The spatial trace of the connection is 

which receives the first and second order perturbations 

r 

and 

a{r)9 {2) = -g l ^ l> g im , T + ^g l3{2 



0(5 2 ) 



= ^ ij om),r = -^v 2 $ (1) 



0(5e 



(6.9) 



(6.10) 



(6.11) 



1 



i(2) 



= $ 

45 



2-2\ 



0(S 2 e 



Using the Fourier transformation (|4.35[) , taking an ensemble average, and using the result that (V 2 $( 2 )) 
Appendix C), we find from Eqs. (|6.4p and (|6.3p 



(6.12) 

(see 



AH 
IT 



Ma (i8 (V2$(1))2 " + 30 ( V2 *WU 



135£/a 



135ffa 



(6.13) 
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which is consistent with the lowest order result of Ref. Q ■ This spatial average is a boundary term, whose ensemble 
average vanishes. 

Although (|6.13|) vanishes, it contains terms with two more powers of k/ Hq than what one would find in the post- 
Newtonian gauge. It is these terms that Refs. argue will lead to a large effect at higher order in perturbation 
theory. In other words, using the synchronous gauge and defining acceleration in terms of spatially averaged expansion 
parameters can lead to a conceivably large correction. This is in contrast to our earlier method, wherein we calculate 
the observable effect, which is very small. Note that our expansion (|A3|) for Dl(z) contains no four-derivative terms 
like those in (|6.13() . 

VII. CONCLUSIONS 

We have computed the inhomogeneity-induced perturbations to the redshifts and luminosity distances that a comov- 
ing observer would measure to first post-Newtonian order, i.e. we have computed z and HqDl to order e 3 ~ (v/c) 3 , 
and to second order in the density perturbation 5 = (p— (p))/(p). Assuming a flat and matter-dominated background 
cosmology, the perturbed luminosity distance- redshift relation is given by Eq. (|4.31|) . The perturbations to Dl(z) 
depend on the correlation between the peculiar velocities at the observer and at the source. Roughly speaking, these 
perturbations are of order AD^/D^ ~ 10~ 5 when z ~ 0.1. The luminosity distance-redshift relation was averaged 
over viewing angles and over an ensemble of realizations of the density perturbation. The result is gauge invariant, as 
it corresponds to a measurable quantity. We then fit this function to what one would expect in a homogeneous FRW 
cosmology, containing dust and a cosmological constant, to deduce the corresponding perturbation to the inferred 
cosmological constant density. 

The inferred fl\ depends on the limiting redshifts z m i n and z max of the sample, and we summarize the best-fit 
values of for different values of these limiting redshifts in Figure [1] and Table Q] These ensemble averaged results 
indicate that we are justified in fitting low-z supernova data to homogeneous models, as long as we use supernova 
data that spans a large enough redshift range. For instance, assuming that we have luminosities and redshifts from 
Zmin = 0.02 out to z max = 0.15, the errors induced by the "fitting problem" are small: Ha ~ 0.004. Such errors are 
not large enough to explain the measured value Oa ~ 0.7. This is what we would expect, since we have other evidence 
to suggest that our universe contains dark energy from large scale structure surveys, from the CMB power spectrum, 
and from weak lensing. 

In contrast to the small value of the best-fit 17a for the ensemble averaged luminosity distance-redshift relation, we 
find that relatively large errors are possible due to fluctuations in Dl(z), specifically from terms that are linear in 
peculiar velocities. This effect was noted in Ref. [4(| and then calculated in Ref. [12| for an actual nearby supernova 
data set. Wc find that the associated variance in Oa has two components, one that depends on the number of sources 
N, a\ ~ (100/A0O ma:r /0.2)- 6 , and one that does not, a\ ~ 0.03(z maa; /0.2)- 8 . 

It should be stressed that our goal in this paper was only to find a rough estimate of the fitting effect. One 
potential weakness of our analysis is that we have assumed that S < 1, and thus we do not address the effects of 
highly nonlinear structures. Such nonlinear modes could be included by using the full nonlinear power spectrum 
from N-body simulations [38[ , and we estimate that this would change the result by approximately a factor of two. 
Furthermore, we have assumed that the observer is in a random location in the Universe, and has no knowledge of 
his/her own peculiar velocity. One can redo the calculation for an observer who knows and corrects for this velocity. 

It has been claimed that there exists a perturbative instability, where successive orders in an expansion in powers 
of S do not get smaller [(| @, Q. We do not see any indications of such an instability with our method. When one 
defines "acceleration" in terms of only directly observable quantities, as we did in Sections II through V, the fitting 
effect one obtains is small. 
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APPENDIX A: COMBINING THE REDSHIFT AND LUMINOSITY DISTANCE RELATIONS 

Adding the redshift equations l[420|) . (|4~28)) , and (|4~29|) yields 



2r r 2 4r 3 , 4 . 



+ 



2 i 6(1)dr ' + §* o(l) 



+ |<(2) - <(2) + *o(2) - *a(3) + |f ("1(2) ~ <(2)) + \ (fs(l) ~ <£(!)) + 

«e(i) fc (i)+^(i) fc ( i)J s + ^*o(2) 



-— * s (2) - —^9i + — + — 



K(i)) 2 - K(d)' 



2r 
W 



o{l) V s(l)"o(l) 



+®o(i)<(i) + *»(!)<(!) + $ o(i)< ( i) - 3$ s( i)< (1) + + O (e 4 <S 2 ) \ + O (e5 3 ) , 



(Al) 



where the right hand side is evaluated at r = r(A) = —A and t = t(A) = to + A. To point out a few of the above 
effects, the terms linear in velocity and linear in $ correspond to the Doppler effect and the gravitational redshift, 
respectively. We also see the second order Doppler shift with the v 2 terms, and the integrated Sachs- Wolfe effect with 
the integrated terms. The perturbed luminosity distance is found from Eqs. (|4. 2 1 [) . (|4.25|) . and (|4.26|) to be 



D L (\6,<f>) 



(1 + z) i 2r 
Hn 3t' 



I * / (r ») 2 V^ (1) *» + 0( e = 5 ) 
Jo r Jo 

jf (.") 2 v 2 <i> (2)( i," + o( £ ^ 



O (e<5 3 ) 



(A2) 



Here we can see the effects of weak gravitational lensing. Note that as the cosmological portion of the redshift goes 
to zero, and hence r — > 0, the luminosity distance also goes to zero, as expected. 

By combining Eqs. (|A1[) and (|A2|) . we can eliminate A and compute as a function of z, 9, and <p. This 
computation can be carried out explicitly by using the fact that the expressions are power series in e and 5. This 
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procedure gives: 




where, to leading order, d/dz w (3to/2)d/dr. The functions of r and t that appear on the right hand side of Eq. 
(|A3j) are evaluated at r — z/Hq and t = to — z/Hq. Note that the redshift z here is the full redshift as measured by 
the observer. Next we need to average Dl(z, 9, </>) over viewing angles in the observer's rest frame, and also take an 
ensemble average. In doing so, the averages of first order quantities will vanish. We also will find that we will only 
need the second order velocities and potentials to Newtonian order, so that we may compute the lowest-order effect. 



APPENDIX B: NEWTONIAN SECOND-ORDER PERTURBATION THEORY 



In terms of comoving coordinates r = x/a(t) [111 ], the equations of Newtonian hydrodynamics are 

85 I 
dt a 



+ V[(l + *)v„] = 0, (Bl) 



+ -v p + - (v p ■ V v p = p - , B2 

at a a a 

and 

V 2 $ p = 4irp a 2 S , (B3) 

where v p = vm + V(2) + ■ • • is the peculiar velocity, <& p = $n) + ^(2) + ■ ■ ■ is the perturbation to the Newtonian 
gravitational potential, the density contrast is 8 — [p(r,t) — Pa{t)]/po{t), and the zeroth order quantities are given 
in Section II. The Newtonian first order results are very well known; for a detailed review, see Peebles For a 

Newtonian analysis to second order in 5, see Ref. 42]. 

The first order result is that the density contrast consists of mode that grows with time, and one that decays with 
time: 

8 (1) {v,t) = f{v)t^+g{v)t^ , (B4) 

where / and g are functions of the spatial coordinates. We will only consider the growing mode. It is useful to rewrite 
the hydrodynamic equations in terms of their Fourier modes. Writing 

' r " k <V' k r (B5) 



(27T) 
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and 



Eq. (|B3[) becomes 



<T> " = / Wf^' r ' (B6) 



fc 2 $ k = 47r / 9 a 2 5 k . (B7) 

The second order density contrast is 

Qf4 / 5 \ 



this result came from perturbing Eqs. (|Bl|l - (jB3|) to second order and then solving these by using the first order 
solutions, Eqs. (|B4[) and (|B6[) . It can be seen that the expected value of 5( 2 ) vanishes by substituting the mode 
expansion of $m into Eq. ()B8[1 : (5(2)) = 0. We also see from Eq. (IB7|) that depends only on boundary 

conditions; we can choose to add overall constants to $ at each order in 6, and it is natural to choose these constants 
to satisfy = ($ (2) ) = 0. 

Assuming that we only have the growing mode solution of Eq. (|B4|) . we find that the first order peculiar velocity 
is related to the Newtonian potential, 

V(i)M) = -4yVS (1) = -i 1 /3 t V3 v$(1) . (B9) 

This averages to zero but its square does not. The second order velocity perturbation is 

3t 3 

v h) = - i4^3*(i),ii*(i)J ( B1 °) 
which also averages to zero: (i>(2)) = 0. Note that these averages are ensemble averages, not spatial averages. 



APPENDIX C: AVERAGING THE LUMINOSITY DISTANCE- REDSHIFT RELATION 



Now we can scrutinize the terms of Eq. (IA3|i . so that we may find their angular and ensemble averages. Note that 
the angular averages will be performed with respect to the observer's angles (9,4>), and so we will need to use the 
Jacobian given in Eq. I|4.13j) . The first three terms of Eq. (IA3|) only depend on the background cosmology, and are 
unchanged after averaging, and all terms that are to first order in S will have a vanishing ensemble average. As shown 
in Appendix B, terms that depend on uLn and $(2) also average to zero. 

In addition, there are many terms that have vanishing ensemble averages because they contain an odd number of 
spatial derivatives of the potential, such as 



< { l)*o(l)) = {< {1) ®s(l)) = , (CI) 

^)|X(i))= > ( C2 ) 

(4)<w} = > ( C3 ) 

^41^')=°' (c4) 

et cetera. We also find that 

+«<«!) kfi)) s - (wfl(i)*fi) k (i)) J ~0(e 4 ) , (C5) 



25 



since ~ ^(1)^(1) ~ e3 ' anc ^ taking the difference of the averages at the source and at the observer introduces 

another factor of z ~ e. 

We can further rewrite the average by exploiting the power law scaling v?^ oc i 2 / 3 , to find 



<<x>> = ((^V/ 3 V$ (1) ) 2 }«((*rV$ (1) ) 2 }(*o-r) 
= <^ (1) )[l-z + 0(z 2 )] . 
We also use ((«L) 2 ) = (v/U)^, and introduce the two point correlation function f(r) 

(<(!)<(!)) = (l-^) + , 

where /(r) is defined by 

(n • v(r , t)n • v(r + rn, t)) = \(v 2 o{1) ) [1 + /(r)] , 



2/3 



(C6) 



(C7) 



(C8) 



and n is a unit vector that defines the viewing direction. 

We can write this correlation function in terms of a more general correlation function Cij(r), using the Fourier 
transform of Eq. (|B6j) and Eqs. ip~3"7|) - ([O0"|) : 



(«o(i) ) c ijW = (Vi(ro,to)vj(r +r,t )) 
This function can be rewritten as 

( u o(i)) c y( r ) = H o 



Ait 



2 roo 



d 3 kk l k j A 2 {k)e- 
W 



ik-r 



3 r z 



where 



and 



A ^ = 8^ 



d 3 kA 2 (k) 
P 



l-(k-rr 



dkA 2 (k) 
P 



b'o(fcr) + j2(A:r)] 







00 d 3 kA 2 (k) 
k 5 



3(k-r) 2 - 1 



dkA 2 (k) 
P 



-J2(At)] 



and where we are using spherical Bessel functions of the first kind: 

sin a; 



3o{x) = 



x 



and 



It follows that 



sm x 5- cos 2: 

ar 



(C9) 



(CIO) 



(Cll) 



(C12) 



(C13) 



(C14) 



^«i))[l + /(r)] - (v 2 o{1) )nVc l3 (r)=H 2 

r<2i.2 poo 

= / ydyT 2 (y) 



-A(r)+B(r) 



1 . / k eq zy 



3 j2 



keqZy 



(C15) 



where k eq = 1/A C ~ 10 3 £fo- We plot l+/(r) in Figure[2j we see that it falls to approximately zero for r ^S> X c ~ 10 Mpc, 
and thus we do not expect it to be important when measuring the distances to supernovae at redshifts z ~ 0.1. Note 
also that / becomes negative for large enough r. 
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Using these simplifications, we finally get 



D L (z) = JL(\ + 1 Z - iz 2 j + AD r L ™(z) + AD c r r (z) , 

where AZ?£ ms (z) is the perturbation that depends on RMS quantities at a given point, which vanishes: 

AD r L ms (z)=0 , 

and ADf rr (z) is the perturbation that depends on /. To sublcading order, this is 



AD c L orr (z) 



(l + *) 2 (<p) 
Ho 



~f(- , 



1 /'(-^](l-2z) 



O 



ffS 2 



(C16) 



(C17) 



(C18) 



where the subleading terms are suppressed by a factor of X c Hq or z. We will only use the lowest order piece, 



AD L (z) = AD c L orr (z) w - 
C 2 k 3 



Ho 



^o 4 7o 



, T 2/ n ycos(k eq zy/H ) sin(k eq zy/Ho) cos(k eq zy/H Q ) 
" J k eq z/H (k eq z/H ) 2 (k eq z/H ) 3 y 

ain(k eq zy / H ) 



+ 6 



(k eq z/H ) 4 y 2 



(C19) 



APPENDIX D: TRANSFORMING FROM THE STANDARD POST-NEWTONIAN GAUGE TO THE 

SYNCHRONOUS GAUGE 



In the standard post-Newtonian gauge discussed in Section II, we can rewrite the metric in terms of conformal 
coordinates, 



ds 2 = a 2 (ry) [- (1 + 2$ (1) + 2$ (2) ) d V 2 + (l - 2* (1) - 2$ (2) ) SydX'dX*] , 



(Dl) 



where we will only need this to Newtonian order, and now the scale factor is 01(77) = (i]/r]o) 2 . We will define 770 « 3to 
to be the conformal time today. This new time coordinate is related to that of Sections II - V by 



77 = 3 



and the radial coordinates are related by 



"-'5 



- 2 / 3 r r 2 / r 4 



-2/3 r 2 / 4 

1+ 9? + °U 



-i + O (te 2 



- + O (re 2 
a 



(D2) 



(D3) 



where J2 = ^(X 1 ) 2 + (X 2 ) 2 + (X 3 ) 2 . Thus, we see that the potentials are the same as before, to Newtonian order, 
except that they now are in terms of comoving distance X % and conformal time 77. We also now use Cartesian 
coordinates for simplicity. 

Our goal is to transform to the synchronous gauge, with new coordinates x M = (r, x l ), where the line element has 
the form 



ds 2 = a 2 {T)g y , v dx ti x v = a 2 (r) [-dr 2 + g io dx l x°] . 
In this gauge, g TT = — 1 and g T i = gi T — 0. We make the following ansatz for the new coordinates: 

V = t + ( T , x) + / (2 ) (t, x) + O (r £ 4 ) 

and 

X i = x 1 + h\ 1} (t, x) + h\ 2) (r, x) + O (Pe 2 ) , 



(D4) 
(D5) 
(D6) 
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where hL*. ~ 8x % , /i? 2 % ~ (5 2 x% /(!) ~ 5tqE 2 , /( 2 ) ~ 8 2 toS 2 , and To ~ 770 is the time today. We are also assuming that 

we have the growing mode only, for which we have the power law scalings oc r° and $(2) oc t 2 . 
In order to find the new metric, we will need the relations 



a 2 (v) = a 2 (r) 



1 + 7/(1) + 7/(2) +0(e 4 ) 



and 



4> (1) (r,, X) + 4> (2) ( V , X) = (r, s) + $ (2) (r, 5) + + 0(e 4 ) + 0(S 3 ) 

Using these and the coordinate transformations (|D5[) and (|D6[) . we find 



implying 



and 



ffrr = - I 1 + -/(I) + -/(2) + 2$ (1) + 2$ (1) ,i/l l (1) + 2$ (2 ) + 2/ (1) + 2/ (2) ) + = -1 



-/(i) + + = o 



-/( 2 ) + 2$ (1) ,^ (1) + 2$ (2) + 2/ (2) - = . 

Similarly, the time-space component of the new metric is 

9n = -/(!),< - /(a),i + ^(i) + ftp) + + ( e3 ) = 

and this implies 

- + ^(i) = 

and 

- /(2),i + ^(2) + ^(l)/(l) = ■ 

Equations ijDlO]) , l[DTT|l . (|D1~3|) and ((Dl4|) are solved by 

r ^4 
/(i) = -o $ (i) + 3 ' 
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(D10) 

(Dll) 
(D12) 
(D13) 
(D14) 

(D15) 



/(2) = -~$(2) + ^ (V1>(1)) 2 + ^ - ~/4*(l),< . 



(D16) 



and 



'(l) 



h' 1 — 



T 2 T 4 T 2 ■ T 2 



(D17) 
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where the arbitrary constants A and B and functions h l (x) and /iq(x) represent residual gauge freedoms associated 
with synchronous coordinates. Setting A and J3 to zero will give us comoving coordinates. We can imagine comoving 
coordinates to be fixed on some spacelike hypersurface from which the worldlines of freely falling particles emanate. If 
we set all of the clocks carried by these particles to the same time on this spacelike hypersurface, then A = B = 0. The 
residual functions h l and h l correspond to simply changing the coordinates on the spacelike hypersurface from which 



worldlines emanate, and we will set h l Q = h z Q = 0. Using this solution for the appropriate coordinate transformation 
we find the spatial part of the new metric to be 



9ij -— Sij 



+h 



4 4 
1 + 7/(D + -/(2) 



2$ 
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2$ 
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2 < £(i),fc/i(i) 



(2)iJ I" "(2)j 



h(2)j,i + ^(l)fe,i^(l)fe,j 



:/(D - 2 *(D 



+ ^-*(i),ifc*(i) )J -fc + O (e 2 ) + 0(^ 3 ) 
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